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Preface 
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i. 
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E 


Inertial  navigation  systems  and  inertial  instruments 
have  improved  so  much  since  their  introduction  that  today  the 
unmodeled  gravity  anomalies  are  becoming  the  dominant  error 
source  in  the  navigation  error  equations.  This  thesis  pro- 
poses the  use  of  precision  range  and  angle  measurements  be- 
tween two  accelerometer  triads  and  an  Extended  Kalman  Filter 
to  estimate  the  anomalous  gravity  gradient. 

Prior  to  undertaking  this  study,  my  knowledge  of  geodesy 
and  nonlinear  estimation  theory  was  minimal.  This  study 
presented  a difficult  challenge,  but  I was  aided  by  the  out- 
standing teaching  abilities  of  Dr.  Peter  S.  Maybeck  and  by 
the  guidance  and  support  of  my  advisor  Captain  J.  Gary  Reid. 

To  these  two  men  and  to  Major  Salvatore  Balsamo,  whose  advice 
and  humor  helped  at  critical  points,  I owe  a great  debt. 

The  true  value  of  this  thesis  is  not  to  be  found  in 
these  pages,  however,  but  in  the  increased  knowledge  and 
analytical  abilities  I take  with  me  as  a result  of  this  effort. 
It  is  this  last  for  which  I am  most  grateful. 


Richard  W.  Smart 
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• Abstract 

A preliminary  feasibility  study  was  done  on  a proposal 
to  estimate  the  anomalous  gravity  gradient  via  precision 
measurements  between  two  accelerometer  triads  moving  through 
a gravity  field.  The  state  equations  to  describe  the  system 
were  developed,  and  an  Extended  Kalman  Filter  was  designed  to 
exploit  the  measurement  data.  A simple,  first-order  numerical 
approximation  was  used  in  gradient  estimation.  In  addition, 
a brief  conceptualization  of  the  hardware  which  might  be  used 
i r\  an  operational  system  of  this  type  is  included.  The  results 
indicate  that  the  numerical  ill-conditioning  in  the  problem 
makes  the  Extended  Kalman  Filter  a poor  choice  of  estimator. 
Even  after  scaling  variables  and  employing  a U-D  Covariance 
factorization  update  algorithm,  numerical  instabilities  in 
the  filter  were  observed.  Recommendations  for  further  work 
on  this  topic  are  included. 
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ESTIMATION  OF  THE  ANOMALOUS  GRAVITY  GRADIENT 
VIA  PRECISION  ANGLE  TRACKING 

I . Introduction 
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Statement  of  the  Problem 

State-of-the-art  in  terrestrial  inertial  navigation  has 
reached  the  point  where  one's  imperfect  knowledge  of  the 
earth's  gravitational  field  (G,  due  to  mass  attraction  alone) 
is  assuming  the  dominant  role  in  the  error  equation  (Ref. 
17*269  and  Ref.  bi2) . In  particular,  INS  errors  due  to 
gravitational  anomalies,  or  to  the  unmodelled  changes  in  the 
gravity  gradient,  are  predominant  (Ref.  17*269). 

This  thesis  addresses  the  feasibility  of  estimating  the 
gravity  gradient  using,  among  other  measurements,  precision 
angle  tracking  between  two  high  accuracy  IMU.  Other  methods 
of  estimating  the  gravity  gradient  such  as  through  the  use  of 
gravity  gradiometers  now  under  development  will  not  be  dis- 
cussed. Should  feasibility  be  shown  in  this  study,  the 
potential  benefits  could  include  more  accurate  inertial  navi- 
gation performance  and  an  alternative  to  gradiometers  for 
rapid  and  wide-area  gravity  surveying. 

Of  the  several  estimation  techniques  ir  use  throughout 
the  scientific  and  engineering  communities,  the  Extended 
Kalman  filter  (EKF)  has  been  chosen  for  this  study.  This 
choice  of  the  EKF  estimator  was  based  on  the  non-linear  dyna- 
mics involved  in  the  problem  and  a desire  to  keep  the  esti- 
mator as  simple  as  possible. 
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Objectives  of  the  Study 

The  primary  objective  of  this  study  is  to  determine  the 
feasibility  of  using  precision  range  and  angle  tracking  be- 
tween two  IMUs  and  employing  an  EKF  to  estimate  the  gravity 
gradient  in  a moving  vehicle  environment.  To  investigate 
feasibility»  the  problem  is  broken  down  into  the  following 
parts t 

1 . Develop  the  hypothetical  design  of 
the  hardware  configuration. 

2.  Develop  a truth  model  and  an  EKF 
estimator. 

3.  Generate  a "flight"  profile 
representing  a typical  mission. 

w 4.  Using  the  models  developed  in  (2) 

and  the  flight  profiles  of  (3)»  carry 
out  a Monte  Carlo  analysis  of  the 
gradient  estimation  technique. 


Assumptions  and  Limitations 

Because  this  is  a feasibility  study,  many  simplifying 
assumptions  will  be  made.  Also,  optimistic  performance  char- 
acteristics of  the  various  measurement  devices  will  be  used. 
These  simple,  optimistic  assumptions  are  justified  in  that 
they  simplify  the  analysis  and  at  the  same  time  they  test  the 
hypothesis  that  this  technique  can,  indeed,  accurately  esti- 
mate the  gravity  gradient.  Failure  of  the  process  to  work 
under  optimistic  assumptions  would  clearly  indicate  that 
further  investigation  along  these  lines  would  not  be  fruitful 
The  first  assumption  made  is  that  the  vehicle  (aircraft, 
ship,  or  submarine)  used  in  the  gradient  survey  will  be  op- 
erated in  a benign  environment  (calm  air  or  sea)  and  it  will 
not  be  undergoing  violent  manuevers.  Nevertheless,  there 
will  be  vibration  and  vehicle  bending  modes  to  contend  with. 

In  this  study  it  is  assumed  that  vibration  and  bending  effects 
can  be  minimized  by  careful  design  and  construction,  and  these 
effects  are  not  considered  further. 

In  the  flight  profiles  tested  only  straight  and  level 
flight  is  simulated.  This  would  typically  be  true  for  a 
mapping  mission  such  as  envisioned  here.  This  assumption  is 
really  not  too  restrictive  since  accelerometers  (here  assumed 

^ Q p 

sensitive  to  10  m/sec  ) can  measure  and  hence  account  for 
vehicle  motion.  Accelerometers  with  thie  sensitivity  have 
been  demonstrated  according  to  Slater  (Ref.  24i198)  for  earth- 
based  instruments.  Instrument  misalignments  are  assumed 
known  apriori  or  can  be  obtained  through  proper  alignment 


and  calibration  techniques. 

The  system  to  be  studied  is  non-linear  and  therefore 
non-linear  estimation  techniques  must  be  used.  The  EKF  used 
here(  as  with  most  methods  of  treating  non-linear  problems, 
assumes  that  a linear  system  model  operating  about  some 
reference  trajectory  is  an  adequate  representation  of  the 
original  non-linear  system.  Implicit  in  this  assumption  is 
the  premise  that  deviations  from  the  reference  trajectory 
are  small  enough  that  linear  methods  for  treating  them  are 
still  valid.  Furthermore,  in  the  EKF  the  model  is  relinear- 
ized about  the  new  estimated  trajectory  at  each  update. 
Further  discussion  on  the  assumptions  made  in  developing  the 
EKF  can  be  found  in  Chapter  III  of  this  report. 

Only  two  flight  profiles  are  simulated.  The  first  pro- 
file simulates  an  aircraft  flying  west  at  100  mi/hr  at  a 
constant  altitude  of  10,000  feet  and  at  a constant  35°  north 
latitude.  The  initial  conditions  for  the  flight  profile  are 
latitude  35°  north,  longitude  239°  east  measured  with  re- 
spect to  the  Greenwich  meridian.  These  initial  conditions 
were  chosen  since  they  correspond  closely  with  the  center 
coordinates  of  the  point-mass  gravity  model  region  used  to 
generate  MtrueN  gravity  data. 

The  second  profile  generated  simulates  a ship  or  sub- 
marine at  the  surface  of  the  ocean  traveling  due  west  at  6 
knot 8 with  the  same  initial  point  as  in  the  first  profile. 
Although  the  initial  point  and,  indeed,  the  whole  profile  is 
over  land  in  reality,  the  identical  path  was  chosen  for  both 
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simulations  for  convenience  and  for  comparative  reasons. 
Intuitively,  the  slower  and  lower  in  altitude  profile  should 
provide  better  gradient  estimation  accuracy  than  the  higher 
speed,  higher  altitude  simulation.  One  fundamental  reason 
for  this  is  the  attenuation  of  the  gradient  with  increasing 
altitude.  Also,  the  slower  speed  allows  more  integration 
time  for  the  anomalous  gravitational  accelerations  to  affect 
position  changes. 

The  actual  gravity  and  gravity  gradient  data  for  the 
truth  model  was  obtained  from  a point-mass  gravity  model 
included  as  part  of  a general  spherical  harmonic  and  point- 
mass  gravity  model  computer  program  called  GTGRV.  This 
program  and  its  associated  data  files  were  obtained  from  the 
Geodynamics  Corporation,  Santa  Barbara,  California.  Pertin- 
ent details  of  the  structure  of  the  program  and  its  use  can 
be  found  in  Reference  20.  All  computer  programs  used  in  this 
study  were  written  in  CDC  Fortran  Extended  and  executed  on  a 
CDC  6600  computer. 

The  point-mass  model  included  seven  point-mass  sets 
contained  within  the  following  geographic  boundaries 1 
latitude  = 12.5°N»  longitude  « 207.5°- 272. 5°E. 

Distribution  of  the  point-masses  within  the  above  boundaries 
varied  from  spacing  of  5°  for  set  1,  where  (°)  represents 
degrees  of  latitude  or  longitude,  to  spacing  of  0.3i25' 
for  set  7,  where  (')  represents  minutes  of  longitude  or 
latitude  (1°  » 60').  A point-mass  model  of  this  detail 
probably  represents  the  "best"  model  that  can  be  obtained 
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from  unclassified  sources. 

A key  assumption  made  in  this  study  is  that 


where 


re,  g.(c) 


9-03  = 9_(c)  . + o(r)  (i-D 

SPHerewmcs  ^ Wt  #W 


total  gravitational  acceleration  as  a 
function  of  position,  £ , 


SFHftrmomcs  = gravitational  acceleration  as  a function 
of  ] C obtained  from  a spherical  harmonic 


model  of  degree  and  order  8,  and 
%($* = gravitational  acceleration  as  a function 

of  £ obtained  from  a point-mass  model. 
The  s&W  ^j^^^term  in  eqn.  (1-1)  is  used  here  to  represent 


the  anomalous,  or  unmodelled,  gravitational  effects. 

The  hardware  configuration  is  assumed  to  include  an 
electro-optical  tracker  capable  of  measuring  range  to  an 
accuracy  of  10“ -’m  and  angular  deviations  of  the  line-of- 
sight  to  an  accuracy  of  l^Urad.  A triad  of  accelerometers 
sensitive  to  10  iq/sec  is  mounted  at  the  tracker's  position, 
and  an  identical  triad  of  accelerometers  is  mounted  on  a 
target  platform.  Although  the  accelerometers  would  not  have 
to  be  part  of  an  IMU,  in  many  applications  several  IMU's 
would  be  available  in  the  aircraft  or  ship  and  could  be  used 
for  this  purpose.  This  assumes,  of  course,  that  the  instru- 
ments in  the  IMU  are  of  the  sensitivity  mentioned  above.  For 
convenience  in  discussion,  the  two  accelerometer  triads  will 
be  called  IMU's.  The  two  IMU's  are  initially  separated  by  a 
known  distance  (10  meters  in  this  study)  and  are  flown 


through  a known  gravity  field. 

Another  of  the  critical  simplifying  assumptions  made  is 
that  all  accelerometer  and  tracker  noise  sources  can  ade- 
quately he  modelled  as  independent  white  Gaussian  noises. 
Justification  for  this  assumption  rests  on  the  premise  that 
very  precise  and  accurately  calibrated  instruments  would  be 
used,  and  hence  most  of  the  systematic  or  bias  errors  would 
be  calibrated  or  compensated  for. 

A system  to  completely  estimate  all  five  of  the  indepen 
dent  elements  of  the  nine  element  gravity  gradient  matrix 
would  require  a minimum  of  three  target  IMUs.  One  possible 
arrangement  for  a full  system  configuration  is  shown  in 
Figure  1. 


Fig.  1 Full  System  Configuration 


The  trackers  in  Fig.  1 are  needed  in  order  to  estimate 


II.  System  Modelling 

Concept  Overview  and  Chapter  Outline 

The  problem,  as  stated  earlier,  is  to  develop  a means 
of  estimating  the  anomalous  gravity  gradient.  Gradiometers, 
a new  class  of  inertial  instruments,  are  presently  undergoing 
extensive  development  and  test  by  a number  of  manufacturers. 
These  instruments  are  specifically  designed  to  measure  the 
elements  of  the  total  gravity  gradient  tensor  in  real-time 
so  that  compensation  for  their  effects  on  navigation  estimates 
can  be  made  (Ref.  4»34).  However,  since  operational  gradio- 
meters are  no  near-term  certainty,  particularly  for  moving 
vehicle  applications,  alternate  methods  of  performing  grad- 
ient mapping  bear  investigation. 

Although  Einstein's  Principle  of  Equivalence  (Ref.  2«2) 
eliminates  the  accelerometer  as  a gravity  sensor,  two  accel- 
erometers in  the  same  dynamic  environment  with  parallel  sen- 
sitive axes  and  separated  by  a small  distance  can  be  used  to 
sense  differential  acceleration  (Ref.  4i36)»  As  long  as  the 
dynamic  environments  of  the  two  accelerometers  are  the  same, 
the  differential  acceleration  between  them  must  be  due  to 
gravitational  variations  (Ref.  4«36).  The  principle  of  using 
differential  accelerations  to  estimate  the  spatial  gravita- 
tional variations  (i.e.,  anomalous  gravity  gradient)  is  the 
basis  for  this  study. 

( Consider  two  platforms  each  supporting  a triad  of  accel- 

erometers and  separated  from  each  other  by  a short  distance. 
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Further,  consider  these  two  platforms  moving  through  a gravity 
field*  One  platform  also  carries  a tracker  for  reasons  which 
will  be  discussed  later.  The  platform,  accelerometers,  and 
tracker  together  will  simply  be  referred  to  as  the  "tracker." 
The  other  platform  is  termed  the  target  platform. 

The  relative  acceleration  between  the  tracker  and  tar- 


get coordinate  frames  can  be  expressed  as 


where, 


aUl  = -fflr  (2-1) 

where,  = relative  acceleration  vector 

expressed  in  the  i- frame, 
fi  J>  = inertial  acceleration  vector  at  the 
target  expressed  in  the  i- frame,  and 
Pi  = inertial  acceleration  vector  at  the 
tracker  expressed  in  the  i-frame. 
Britting  (Ref.  2 02)  gives  as  an  equation  for  the  output  of 
an  ideal  accelerometer, 

_ OL  ••  I m r-V.  A — 1 


where i 


K' 

K/  * summation  over  all  the  bodies  of  the 
universe  except  the  earth, 

G*  * gravitational  acceleration  at  the 

instrument  location  due  to  the  earth, 

••  | 

£ * inertially  referenced  acceleration, 

* coordinate  transformation  matrix  relat- 
ing the  inertial  axes  i to  accelerometer 


axes  a,  and 


GK'»  gravitational  acceleration  at  the  earth's 


mmmx: 


E I 

C ) 


center  of  mass  due  to  all  of  the  k 
bodies  in  the  universe  except  the  earth. 
Rearranging  eqn.  (2-2)  and  expressing  all  quantities  in 
the  i- frame  yields 

£'  = £!+£5-2Rfc  - Gt.]  (2-3) 

K' 

Equation  (2-3)  actually  represents  two  expressions,  one  for 
ft  Xb  and  one  for  pj  Tt  • Substituting  these  two  equations 
into  eqn.  (2-1)  gives 

and  grouping  like  terms  yields, 

Qs'fl  = (fb~£r)  +(§b~  §r) 

Each  component  of  the  last  term  in  eqn.  (2-5)  represents  the 
difference  between  the  gravitational  effect  at  the  center  of 
the  earth  due  to  all  other  bodies  in  the  universe  and  the 
same  effect  at  the  instrument  location  (either  at  the  target 
or  the  tracker).  ' According  to  Britting  (Ref.  2<32),  these 
difference  terms  for  the  moon  and  the  sun,  which  have  the 
largest  effect,  are  on  the  order  of  10”^|Gl . Heuristically, 
since  these  difference  terms  would  be  nearly  identical  for 
two  bodies  only  10m  apart  (as  assumed  here)  near  the  earth's 
surface,  the  last  term  in  eqn.  (2-5)  becomes  vanishingly 
small  and  can  be  neglected.  This  assumption  should  be 
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(2-4) 


(2-5) 


fi'ril  = (fb  ”fr)  +(&b  “ &t) 


verified  more  rigorously  in  any  subsequent  study  of  this 
proposal.  Equation  (2-5)  can  then  be  rewritten  as 

(2-6) 

The  first  term  on  the  right-hand  side  of  eqn.  (2-6)  accounts 
for  all  the  non- field  forces  acting  on  the  target  and  tracker 
and  the  second  term  accounts  for  the  field  forces,  here 
assumed  due  only  to  the  earth's  gravitational  field.  The 
gravitational  acceleration  vectors,  and  Qt  , can  be 
written  as  the  sum  of  a modelled  gravitational  acceleration 
and  an  unmodelled  gravitational  acceleration,  thus, 

■'  _ •. 


/*•  — • 

T ^ (uwnwdellei) 

Substituting  eqn.  (2-7)  into  eqn.  (2-6)  gives 


unm»<UJIe^ 


(2-7a) 

(2-7b) 

(2-8) 


The  last  term  in  eqn.  (2-8)  can  be  replaced  by  £g.  and 
an  error  term  added  to  account  for  accelerometer  errors. 
Thus,  eqn.  (2-8)  can  be  rewritten  as 

(2-9) 


'mm  • H * - 


a-'*t = (fi-fr)  + 

where  £ represents  the  difference  between  the  accelerometer 
error  vectors  (gb-er)  • 

If.  (1)  fb  and  can  be  measured  in  coordinate 
frames  which  are  of  known  orientation  with 
respect  to  inertial  space,  and 
(2)  Q’fm.t  can  be  estimated  by  an  independent 
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source  of  information* 

then  eqn.  (2-1)  may  be  used  to  estimate,  Sg  , the  difference 
in  unmodelled  gravity  between  the  two  locations.  The  esti- 
mate, in  turn,  can  be  used  to  estimate  the  elements  of  the 
gravity  gradient  tensor  via 


fid  » \f 

1 = [il  ||  ill  (2-l0) 

j \Jkc*  SCy  irmJ 

It  is  critical  to  satisfy  conditions  (1)  and  (2)  above 
for  this  proposal  to  work  satisfactorily.  This  thesis  add- 
resses  only  the  problem  of  estimating  4*1  * In  particular, 
this  thesis  examines  the  capability  of  a line-of-sight  (LOS) 
tracker  producing  range  and  angular  deviation  measurements 
to  adequately  estimate  Qrral  . These  measurements  are  then 
employed  by  an  EKF  to  solve  for  the  % components . 

Full  examination  of  the  first  condition  (maintaining 
precision  orientation  information  about  the  two  accelerome- 
ter triads)  is  beyond  the  scope  of  this  the3isi  however,  some 
potential  approaches  to  this  problem  are  suggested  in  the 
discussion  of  hardware  configurations. 

In  the  remainder  of  this  chapter,  the  foregoing  estima- 
tion problem  will  be  developed  in  detail  and  some  potential 
hardware  configurations  are  described.  Specifically,  the 
next  topic  is  a discussion  of  coordinate  frames,  followed  by 
the  development  of  the  state  equations  for  the  estimator. 
Finally,  the  chapter  is  concluded  by  a brief  description  of 
some  potential  hardware  components  which  might  be  employed 
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in  an  operational  system  as  envisioned  here 
Coordinate  Frames 


Table  I is  a list  of  the  coordinate  frames  used  in  this 


study 


Table  I 


List  of  Coordinate  Frames 


The  Earth-Centered-Inertial  (ECI)  or  i- frame  has  its 
origin  at  the  center  of  the  earth.  The  x^-axis  points  in  the 
direction  of  the  vernal  equinox,  or  the  constellation  of 
Aries.  The  y^-axis  lies  in  the  equatorial  plane  with  the 
x^-axis,  but  90°  east  of  the  x^-axis.  The  z^-axis  points  in 
the  direction  of  the  north  pole.  The  i- frame  is  non- rotat- 
ing with  respect  to  the  stars  and  the  earth  turns  relative 


The  Earth-Centered  Earth- Fixed  (ECEF)  or  e- frame  also 
has  its  origin  at  the  earth's  center  and  the  xft-ye  plane  also 
is  coplanar  with  the  equatorial  plane.  The  z -axis  points  in 
the  direction  of  the  north  pole  and  the  x-axis  lies  in  the 


Frame 

Sub/superscript 

Components 

Inertial  (ECI) 

i 

yit 

Earth  (ECEF) 

e 

V V ze 

Tracker 

T 

Xip,  yT. 

Target  Platform 

b 

V yb’  zb 

Line-of-sight 

LOS 

^OS*  yL0S*  ^OS 

Accelerometer 

a 

V ya’  za 

Greenwich  meridian.  The  y -axis  is  90  east  of  the  x-axis. 
At  time  t=0,  or  the  initial  time  of  the  problem  under  study, 
the  xe-axi8  is  assumed  aligned  with  the  x^-axis.  The  e-frame 
has  axes  fixed  in  the  earth  and  therefore  rotates  about  the 
se-axi8  at  the  earth’s  rotation  rate,  CO  ie=?. 2921151^7  x 10"^ 
rad/sec.  The  angle  between  the  x^  and  the  x@  axes,  then,  is 
found  by  the  product  COiet,  where  t is  the  elapsed  time  since 
the  problem  began.  The  e- frame  is  important  in  this  study 
in  that  the  computer  program  used  to  calculate  the  true  grav- 
itational accelerations  and  gradients  produces  its  values  in 
e-frame  components.  Figure  2 shows  the  relationship  between 


inert" ial  __ 
reference 
meridian 


Greenwich 

Meridian 


Fig.  2 ECI  and  ECEP  Coordinate  Frames 


The  tracker  or  T- frame  has  its  origin  at  the  tracker's 
location  with  the  x^-axis  oriented  along  the  boresight  of  the 


tracker.  The  y,j,  and  zT  axes  then  form  an  orthogonal  triad 
with  x,p.  The  T- frame  is  instrumented  with  three  rate  gyros 
and  three  accelerometers.  The  measurements  are  made  in  the 
tracker's  coordinate  frame  ( T- frame ) . Because  perfect  track- 
ing is  not  possible  in  practice,  the  T- frame  is  not  perfectly 
aligned  with  the  true  line-of-sight  (LOS)  coordinate  frame 
in  which  the  target  lies.  These  two  coordinate  frames  can 
be  related  by  two  Euler  angle  rotations  and  Se  as  shown 
in  Figure  3* 


Tracker 


Fig.  3 Tracker  and  LOS  Coordinate  Frames 


As  can  be  seen  in  Figure  3»  rotates  the  LOS- frame 
about  the  z^-axis,  and  $6  rotates  the  LOS- frame  about  the 
y^-axis.  The  transformation  maxrix  which  will  map  a vector 
in  the  T- frame  into  a vector  in  the  LOS- frame  is  denoted  by 
the  symbol  • For  Si\  and  io.  in  the  microradian  range 

as  would  be  the  case  for  near-perfect  tracking,  Cx*  i® 


4 


given  by, 


Se 


8yv  -fe~ 

\ o 

O I 


(2-11) 


The  target  platform  or  b- frame  is  an  orthogonal  triad 
with  origin  at  the  target's  location.  The  b- frame  is  instru- 
mented with  gyros  and  accelerometers  and  its  attitude  with 
respect  to  inertial  space  is  known. 

The  accelerometer  or  a- frame  actually  represents  two 
frames,  one  at  the  tracker's  location  (a,p- frame)  and  one  at 
the  target’s  location  (a^- frame).  This  frame  represents  the 
three  sensitive  axes  along  which  the  accelerometers  make 
measurements,  since  it  is  assumed  that  there  is  no  misalign- 
ment between  the  instruments  and  their  mountings.  In  this 
application  the  a^  and  b- frames  are  considered  perfectly 
aligned  as  are  the  a^  and  T- frames.  Later  in  this  section 
reference  is  made  to  the  a- frame,  and  it  is  to  be  understood 
to  apply  equally  to  the  a^  and  a^,- frames. 

The  relationships  between  the  i,  e,  and  Tracker  or  T 
coordinate  frames  will  now  be  developed.  This  development 
and  much  of  what  follows  in  this  chapter  closely  parallels 
the  similar  development  by  Mitchell  in  Reference  18. 

From  the  gyros  on  the  tracker,  one  can  obtain  the 
tracker's  attitude  with  respect  to  the  i- frame.  As  shown  in 
Figure  4 (Ref.  18il2),  the  T- frame  can  be  related  to  the  i- 
frame  by  two  Euler  angle  rotations,  here  called  0 and  £$  . 


I 
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As  given  by  Mitchell  (Ref.  I81I3),  the  coordinate  transforma 
tion  matrix  from  the  i- frame  to  the  T- frame  is 

^ ”co*  8 cos  <}6  sin  0 cos  5*5  -s\r\  d> 

C;  = Cos  8 O (2-12) 

cos  8 sin  ^ S»»v8  sinqi  coss^ 


Pig.  4 Inertial  and  Tracker  Coordinate  Frames 

It  remains  to  determine  the  Euler  angles  and 
Again,  borrowing  from  Mitchell  (Ref.  18s 14), 

0 = Tan'(B&^j  Rx>0 

-I^V-SWRx^ 


* = t<xn’’  [(?C£+v^]  ‘ 

where  Rx  • Ry  , and  are  the  three  components  of  the 


relative  position  vector  of  the  target  from  the  tracker  ex- 
pressed in  the  i- frame. 


State  Equations 


Figure  5 is  a free-body  diagram  of  the  gradient  mapping 
system  and  is  used  to  aid  the  derivations  of  the  state  equa- 
tions* In  Fig.  5»  JR  is  the  relative  position  vector  between 
O and  b » the  origins  of  the  tracker  and  target  platform 
frames  respectively. 


Fig.  5 Free-Body  Diagram  of  Basic  System  Elements 


The  magnitude  of  R , or  | P>\  , is  the  scalar  range  vari- 
able r . Taking  the  second  derivative  of  with  respect  to 
inertial  space  and  applying  the  Theorem  of  Coriolis  as  needed 
gives  (Ref.  I8il9)» 

pfB  = p1,£  + 2»*x  x£ 

x m x/Vi>  xR\  (2-16) 


time  derivative 


second  time  derivative , 

time  derivative  with  respect  to  the  i and  LOS- frames 
respectively, 

angular  velocity  vector  of  the  LOS  with  respect  to 

the  i- frame,  and 

vector  cross  product  operation. 


Equation  (2-16)  is  a physical  vector  equation  (i.e.,  the 
vector  components  of  R are  not  expressed  in  any  particular 
coordinate  frame).  To  make  equation  (2-16)  a mathematical 
vector  equation  choose,  say,  the  LOS  coordinate  frame  to 
work  in.  Defining  the  term  p-* as  <Xra\  (the  relative 
acceleration  between  points  O and  b of  Pig.  5)  and  writing 
out  the  vectors  in  eqn.  (2-1 6)  in  component  form  yields. 


rr*rW] 


Solving  eqn.  (2-17)  for  the  derivative  terms  p^r  • 
p tA*.  , and  p . yields  the  following  three  scalar 


equations i 


Puss^wosy  - '"•p  a^a  “ ^^u»x^u>»e  (2-19) 


Puks^ujfta  = ^r^ni y - ^u»sy  ■+'k)u>sxa\Asy  ( 2“2°) 

Physical  insight  allows  one  to  write  a fourth  scalar  differ- 
ential equation  as, 

Pu*r=Vr  (2-21) 

where  in  eqns.  (2-18)  through  (2-21) i 
T~  = range 
Vr  = range  rate 

Qftl  x = (acceleration  of  the  target  along  the  axis 

minus  acceleration  of  the  tracker  along  the 
axis)»  similarly  for  OlnVy  andCXrelt  • 

As  pointed  out  by  Mitchell  (Ref.  18s 21),  the  angular 
velocities  in  the  above  state  equations  will  be  measured  in 
the  T- frame  since  there  is  no  physical  means  of  measuring  in 
the  LOS- frame.  Additionally,  the  angular  velocity  about  the 
line-of-sight,  {OlASx  , has  no  significance  and  must  be  elim- 
inated from  eqns.  (2-19)  and  (2>20).  It  can  be  shown  that 
&)kPSx  can  be  expressed  as  (Ref.  18i21), 

£Ovp«x  “ ^Tx  +* ^l^y  “ Sfi 60-rfc  (2-22) 

where,  6)tx  * ^Ty  • andC^reare  the  inertial  angular  acceler- 
ations expressed  in  the  T-frame.  Substituting,  eqn.  (2-22) 
into  eqns.  (2-19)  and  (2-20)  yields, 

fW°tf>Sy  = " ^*JDSy+  ^UJS »(^x+^7^Ty~^e 

and, 
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fias^uSfe5  "ir^y  “ &u>sy  (^y  +S|ltt¥y  “*  Sc.  ^>Te.)  (2-24) 

The  relative  acceleration  vector,  • is  expressed 

in  the  LOS-frame j however,  the  accelerometer  outputs  for  the 
tracker  are  available  directly  in  T- frame  coordinates  while 
those  of  the  target  are  expressed  in  b- frame  coordinates. 

As  a first  step,  the  transformation  from  target  accelerations 


to  the  LOS-frame  can  be  accomplished  by 


= Ql  a*  (2‘25) 

where  the  transformation  matrix  c£  is  assumed  known.  Tech- 
niques for  obtaining  information  about  Ql  are  discussed  in 
the  next  section.  Now  the  relative  acceleration  vector  Q^. 


can  be  defined  as 


QTr  = 9-1  - a; 


(2-26) 


-IPS 

Finally,  using  the  transformation  matrix  QT  defined  in 

CPS 

eqn.  (2-11)  yields  the  following  expression  for 

fi£  = £rQTr  . (2-27) 


Expanding  eqn.  (2-27),  expressing  it  in  component  form,  and 
substituting  into  eqns.  (2-23),  (2-24),  and  (2-18)  gives 


fipsYr  = ttrx  ~ Sfe 

+ |ji£arx  (Sv^y-****^ 

+ Dr1*  arx  “^uwy  (jWhy ***«.)] 


(2-28) 


(2-29) 


(2-30) 
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Li  :<r.v' . 


(2-31) 


One  additional  modification  to  the  above  three  state 
equations  comes  from  substituting  eqn.  (2-9)  into  eqns.  (2-28)* 
(2-29) » and  (2-30).  This  gives, 

Pu»Vr  = tfVfrx)  +(9k,f 

<9^-30^,  + S3*] 

~2.)/r  Ut^sy  +-“^*^10^  +£e 

+t'¥-S3x  + ^LO*a  (Su^Ty-tea**)] 


(2-32) 


Putt  USg-  -p- [ffby  "fty)  +"(9by + $3)^ 


^Sy^x  + £y 

+£“  $3x  -to  VflSy  (^noHy  " Sa 


(2-33) 


It  is  assumed  here  that  the  tracker  and  target  are  both 
in  the  identical  dynamic  environment  (i.e.,  vibration,  bend- 
ing, and  other  forces  are  the  same  on  both  bodies).  Equa- 
tions (2-31) » (2-32),  and  (2-33)  can  then  be  rewritten  as 

Pwos^r  = ^3  x + r6°iPSy  ‘*"d*\ps*)‘h%%y  “ + £*  ( 2-34) 


60  = - 


^U»Sy  +*  £* 

+ jj^*  $3x  + 

l33L  - Z>fc  6)^  -60^ y corx  + £y 
+L‘^r"  y (Snft>Ty  “ S£ 


(2-35) 


(2-36) 


• • • 


• > ^ •* 


W" 


• ‘ 11  - — I 
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This  assumption  simplifies  the  dynamic  simulation  employed  in 
the  study,  but  it  is  not  overly  restrictive  since  non-field 
forces  are  sensed  by  the  accelerometers  and  could  thus  be 
accounted  for  in  actual  practice. 

The  error  angles,  and  £& , appear  in  the  above  state 
equations  and  are  a result  of  imperfect  tracking.  To  esti- 
mate these  two  variables  with  the  EKF  it  is  necessary  to 
develop  state  equations  for  them.  In  practice,  the  tracker 
will  provide  measurement  information  about  and  $£  , and 
the  necessary  state  equations  for  these  variables  are  given 
(with  approximations)  by  Mitchell  as  (Ref.  18«19) 

Ptos^  = ox - Se coTx  (2-37) 

= ^U>S  y ~ ^Ty  (2-38) 

Six  state  equations  have  been  developed  in  the  preceding 
discussion,  eqn.  (2-21)  and  eqns.  (2-3*0  through  (2-38). 
However,  in  keeping  with  the  objectives  of  this  initial  fea- 
sibility study,  five  of  these  six  equations  were  simplified 
for  application  both  in  the  truth  model  and  in  the  filter. 

In  eqns.  (2-3*0  through  (2-38),  the  first  simplification 
made  was  to  neglect  terms  multiplied  by  &TI  or  For  the 

high  accuracy  tracking  assumed  here  both  and  $€  will 
have  magnitudes  on  the  order  of  10“^  rad  or  smaller.  For  the 
profiles  under  test,  the  terms  with  brackets  in  eqns.  (2-35) 
and  (2-36)  have  magnitudes  near  10-1^  which  are  approximately 
5 orders  magnitude  smaller  than  the  magnitudes  of  y 

and  p^O)^  . Also,  in  eqn.  (2-34)  the  product  terms 
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involving  Sir^  and  fe  have  magnitudes  near  lO"1^  while  fi 

o 

has  a magnitude  near  10“  . With  these  simplifications  plus, 
the  elimination  of  the  ^Tx  terms  (since  torsional  bending 
is  negligible)  yields  the  final  six  equations  of  the  EKF. 


Puosr=^ 

Pus^psy  = -Sgr  -ZVc.cOyBg 
f f 7 

Pu«6>u>»*=:  -~2-Vr  W3 

r r * 

fv^s^V.  — +’V^4' 

P\0S^e  = tOlPSv  “"^Tv  ■'*^5 


(2-39) 


(2-40) 


(2-41) 


(2-42) 


(2-43) 


(2-44) 


where  Wt  through  W 5 are  zero-mean  white  Gaussian  noises 
added  to  account  for  the  neglected  terms.  All  that  remains 
is  to  develop  the  stochastic  process  models  for  the  % 
components . 

In  this  study,  the  estimation  problem  is  formulated  in 
state-space  with  a dynamics  model  satisfying  a stochastic 
differential  equation.  It  is  appropriate,  therefore,  to 
seek  a stochastic  model  for  the  components . Typically 
the  anomalous  gravity  field  or  "disturbance"  model  is  sepa- 
rated into  north  and  east  deflections  of  the  vertical  and  a 
gravity  anomaly  component  (see,  for  example.  Ref.  261 135-1^2) . 
As  shown  in  Fig.  6,  the  vertical  deflections  are  defined  as 
the  angular  deviations  between  the  true  gravity  vector  ^ , 

which  is  normal  to  the  geoid,  and  the  reference  gravity 


- 'V;.  r~ 


vector,  , which  is  normal  to  the  reference  ellipsoid. 

The  geoid  represents  a surface  of  constant  gravitational 
potential,  and  the  difference  in  height  between  the  geoid 
and  the  reference  ellipsoid  surface  is  called  the  geoid  height 
or  undulation,  N. 


Vertical 

de-flection 


Pig.  6 Gravity  Disturbance  Quantities 

(Ref.  4.111) 


Several  stochastic  models  for  the  disturbance  components 
can  be  found  in  the  literature.  Musick  in  Reference  19  pre- 
sents synopses  of  seven  papers  describing  different  statis- 
tical models.  The  most  basic  disturbance  model  assumes  that 
a first-order  Markov  process  adequately  characterizes  the 
disturbance  components.  This  model  has  some  physical  appeal 
as  it  does  account  for  the  spatial  correlation  of  gravity 
disturbances  that  is  observed  in  nature  (e.g.,  see  Ref.  26. 
135-138).  Carlson,  as  reported  in  Reference  26.136,  has 


suggested  that  an  appropriate  model  is  the  sum  of  a first- 
order  Markov  process  plus  a random  bias.  Other  models,  such 
as  the  second-order  Markov  model  proposed  by  Kasper  (Ref. 

13) 1 are  higher-order  models  which  attempt  to  model  more 
closely  the  empirical  disturbance  data.  The  disadvantage, 
of  course,  in  using  the  higher-order  models  in  an  on-line 
estimator  is  the  increase  in  state  vector  dimension  as  the 
order  of  the  disturbance  model  increases.  For  this  reason, 
a first-order  Markov  process  model  was  chosen  for  use  in  this 
study,  although  any  subsequent  work  to  this  study  should 
investigate  the  benefits  to  be  derived  from  the  higher-order 
models.  It  may  also  be  quite  appropriate  to  use  different 
statistical  models  for  the  downrange,  crossrange,  and  vert- 
ical directions. 

The  first-order  Markov  process  can  be  generated  by  pass- 
ing white  noise  through  a first-order  lag  as  shown  in  Fig. 

7*  The  autocorrelation,  ¥xx.  and  power  spectral  density, 
functions  for  this  exponentially  time  correlated  pro- 
cess are  shown  in  Fig.  8. 

In  reality  gravity  disturbances  are  spatially  correlated 
and  not  intrinsically  time  correlated.  Nevertheless,  one  can 
transform  a spatially  correlated  process  into  a time  corre- 
lated process  and  vice  versa  by  relating  time  and  distance 
through  the  vehicle's  velocity  variable  (assuming  constant 
velocity)  (i.e.,dBvt*  » where  d is  distance,  V is  velocity, 
and  “t*  is  time).  If  this  is  done,  then  the  expression  for  the 
autocorrelation  function  shown  in  Fig.  8 becomes, 


Fig.  7 First-Order  Markov  Shaping  Filter 

(shaded  box  convention  is  used  in  Ref.  6) 


Aoto  correlation 


Power  Spectral  Density 


Z<rx/T 

af  + Q/rf 


Fig.  8 Autocorrelation  and  Power  Spectral  Density 
Functions  for  First-Order  Markov  Process 
(Ref.  13iChap  4,  p.  85) 


where  D is  the  correlation  distance  and  <Tx  is  the  standard 
deviation.  Widnall  and  Grundy  (Ref.  26 i 137)  give  the  follow- 
ing values  for  the  standard  deviation  and  correlation  distance 
for  each  of  the  three  disturbance  components  for  the  western 
U.S.  assuming  zero  correlation  between  components. 


Table  II.  Gravity  Variations  for  Western  U.S. 


Component 

D ( nmi . ) 

east-west  deflection 

26 

10 

north-south  deflection 

17 

10 

anomaly 

35 

60 

(magnitude) 

Up  until  this  point,  all  of  the  gravity  disturbance 


models  discussed  are  meant  to  model  the  disturbance  along 
one  particular  trajectory.  The  situation  in  this  problem, 
however,  is  to  model  the  difference  in  disturbance  between 
two  trajectories  only  10  meters  apart.  Since  differential 
disturbance  models  are  not  available  in  the  literature,  it 
was  decided  to  do  a spectral  analysis  of  the  empirical 
gravity  data  and  to  propose  an  appropriate  shaping  filter  to 
generate  a process  with  the  same  statistics. 

The  power  spectral  density  (PSD)  plots  generated  from 
the  limited  data  samples  available  can  be  found  in  Appendix 
C along  with  a discussion  of  how  they  were  generated.  As  a 
result  of  this  analysis,  a first-order  Markov  process  model 
was  chosen  for  the  10,000  foot  altitude  case. 

The  state  variable  equations  for  the  three  disturbance 
components  based  on  a first-order  Markov  model  can  now  be 
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written  as, 


+ \di(+) 

(2-46a) 

+ Wr(f) 

(2-46b) 

+ w*(t) 

(2-46c) 

where  i *Tx  »*Ty  , andTg  = correlation  times  for  the 

S$x  . &jy  t and  processes 
respectively,  and 
= zero-mean  mutually  independent 
white  Gaussian  noises  with 
statistics 

E £ w,  (+)  vjJ  (T+r)}  = (t)  s (y) 

E £wz(+)  wl  (t+r)}  s cp(A  S(t) 

E tf)  Stf) 

Figure  9 is  a plot  of  the  true  components  for  the 

10,000  foot  altitude  simulation.  The  curve  was  plotted  from 

data  obtained  from  GTGRV.  It  represents  the  gravitational 

accelerations  due  only  to  point  masses  in  the  GTGRV  model 

for  a vehicle  at  10,000  feet  flying  at  a speed  of  100  mi/hr. 

The  data  was  obtained  at  10  sec  intervals.  The  sawtooth 

effect  in  the  curves  is  a result  of  the  sampling  scheme 

employed.  Note  that  the  magnitudes  of  the  So.  components  for 

•I 

the  time  period  of  interest  are  <2.XlOfig*l,  considerably 
smaller  than  the  disturbance  values  given  in  Table  II.  The 
reason  for  this  is  that  the  data  in  Fig.  9 represents  the 
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differential  disturbances  for  two  objects  only  10  meters 
apart,  whereas,  Table  II  values  are  for  one  object.  These 
small  disturbance  magnitudes  contributed  to  the  problems 
which  attended  this  study. 

As  plotted,  the  Sg.  components  are  expressed  in  L0S- 
frame  coordinates  al though  they  were  initially  output  by 
GTGRV  in  e-frame  coordinates.  The  coordinate  transformation 
matrix  used  in  this  case  is 


(2-47) 


LOS 


Although  in  reality  Cg,  is  not  time-variant,  for  the 
short  duration  of  this  simulation  (745  sec  for  the  10,000 
foot  altitude  case)  compared  to  the  dynamics  of  the  system, 
eqn.  (2-47)  remains  valid.  It  is  hoped  that  a first-order 
Markov  model  is  an  adequate  model  for  the  Sg  components  in 
Pig.  9. 

In  contrast  to  the  obvious  periodicity  of  the  Sg.  compon- 
ents in  the  high  altitude,  high  speed  (relative  to  6 knots) 
case  of  Fig.  9»  the  Sg.  components  for  the  0 ft.  altitude,  6 
knot  case  are  plotted  in  Fig.  10.  Note  that,  although  one 
knows  in  reality  that  Sig.  in  this  case  is  also  time- correlated, 
the  duration  of  the  simulation  is  much  less  than  the  correla- 
tion time.  In  view  of  this,  it  was  decided  to  change  the 
model  for  the  % components  in  the  0 ft.  altitude  case  from 
a first-order  Markov  model  to  a random  bias  model,  although 
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the  original  model  is  still  valid . A random  bias  process  can 
be  generated  as  the  output  of  an  integrator  with  no  input 
other  than  an  initial  condition.  The  describing  differential 
equation  is  given  by 

x(+)  =:  0 (2-48) 

In  the  actual  simulation,  a pseudonoise,  a low  strength 
white  Gaussian  noise,  was  added  to  this  model  to  insure  that 
the  EKF  would  track  the  slowly  changing  bias. 

The  shaping  filter  for  this  process  is  shown  in  Fig.  11 
and  the  autocorrelation  and  power  spectral  density  functions 
for  this  process  are  shown  in  Fig.  12. 


Fig.  11  Shaping  Filter  for  Random  Bias 
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Pig.  12  Autocorrelation  and  Power  Spectral 
Density  Functions  for  Random  Bias 
(Ref.  15*  Chap  4,  p.  85) 

O 

Measurement  Equations 

Only  three  measurements  are  assumed  in  this  study  - 
measurements  of  range  and  the  tracking  angular  deviations* 
and  £<l.  In  keeping  with  the  simplifying  assumptions 
made  throughout  this  study,  the  range  measurement  is  assumed 
to  he  the  true  range  value  corrupted  by  a white  Gaussian 
noise.  Thus, 

r«-rTf0t+v,  (2-49) 

where,  is  a zero-mean  white  Gaussian  noise  with  standard 
deviation  of  lxlO'^m.  A more  realistic  measurement  model 
would  depend  upon  the  type  of  instrument  used,  but  for  the 
purpose  of  this  generic  study  eqn.  (2-49)  suffices. 

The  measurements  of  Sn.  and  Sc.  will  not  truly  be  mea- 
surements at  all  in  the  sense  that  they  represent  a true  value 
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plus  some  corruption.  It  will  be  remembered  that  Stv  and  Sc 
were  equal  to  zero  for  all  time  in  the  truth  model.  Instead, 
$i\  and  Sc  "measurements " are  generated  as  samples  from  white 
Gaussian  noise  processes,  each  with  a standard  deviation  (1«“) 
of  lyurad.  Thus, 

Sn  - V2.  (2- 50a) 

Se  = V3  (2- 50b) 

The  values  for  system  parameters,  °*Ty  and  • were 

_7 

generated  by  adding  white  Gaussian  noise  with»  =10  to  the 
true  values  of  and  m . 

The  assumption  here,  of  course,  is  that  the  angular  rate 
sensors  on  the  tracker  produce  measurements  that  can  be 
expressed  as, 

0)Ty  = 0*v>sY  + V+  (2-51a) 

(a*tb  - 0l>  4*  VS  (2- 51b) 

Once  again,  a very  simplistic  model  for  angular  rate 
errors  is  assumed.  It  should  be  recognized  that  more  detail- 
ed modelling  would  be  necessary  in  an  evaluation  of  a partic- 
ular hardware  configuration. 

Hardware  Configurations 

This  discussion  will  include  a general  description  of 
some  of  the  system  components  and  will  illustrate  some  poten- 
tial measurement  techniques.  No  attempt  is  made  to  present 
a detailed  hardware  design,  however,  as  that  is  beyond  the 
scope  of  this  study. 

To  begin  with,  the  range  measurement  is  important  to  the 
overall  accuracy  of  the  gradient  estimation  process  because 
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of  its  appearance  in  eqn.  (2-10).  One  possible  source  of 
relative  range  measurements  is  an  interferometer  shown  in 
Pig.  13. 


Target 

mirror 


TaTionary 

mirror 


Laser 


In  Fig.  13.  a laser  produces  a narrow  beam  of  monochro- 
matic light.  At  the  beam-splitter,  one-half  of  the  incident 
light  is  reflected  toward  the  target  mirror  and  one-half 
toward  the  stationary  mirror.  If  the  two  mirrors  are  each 
an  integral  number  of  wavelengths  from  the  beam-splitter,  the 
reflected  light  arrives  at  the  detector  in  phase  and  a bright 
spot  is  detected.  If  on  the  other  hand,  the  target  mirror 
moves  one-half  of  a wavelength,  then  destructive  interference 
occurs  and  no  light,  or  very  little  light,  reaches  the  detec- 
tor. Thus,  the  intensity  variations  of  light  at  the  detector 
can  be  used  to  infer  relative  movement  between  the  target  and 
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the  beam-splitter.  The  device  as  just  described  would  be 
very  sensitive  to  changes  in  the  characteristics  of  the  op- 
tical path  medium  (e.g.,  changes  in  the  index  of  refraction), 
thermal  or  other  distortions  in  the  mirrors  and  beam-splitter, 
and  vibration.  Evacuating  the  light  pathways  would  alleviate 
problems  of  the  first  kind.  Proper  environmental  controls 
and  manufacturing  processes  could  minimize  problems  of  the 
second  kind.  However,  isolating  the  device  from  vibration, 
especially  in  a moving  vehicle  such  as  an  aircraft  or  sub- 
marine, represents  a formidable  challenge  and  may  be  the 
limiting  factor  in  an  operational  system. 

It  remains  for  the  tracker  to  measure  the  misalignment 
angles  between  the  LOS  and  T-coordinate  frames.  This  could, 
perhaps,  be  performed  by  the  same  instrument  designed  to 
obtain  the  relative  range  measurements  discussed  above.'  One 
interesting  technique  for  accomplishing  this  is  described  in 
Reference  8.  The  instrument,  called  STARTRAC,  provides 
angular  deviation  measurements  between  star  LOS  and  stable 
member  mirror  normal  in  two  degrees  of  rotational  freedom 
(Ref.  8iSection  1,  p.4).  A simplified  diagram  of  STARTRAC 
is  shown  in  Fig.  14.  f 

In  Fig.  14,  the  stable  member  mirror  could  be  mounted 
on  the  tracker  with  the  normal  to  its  surface  being  colinear 
with  the  Xq,-axis.  A light  source  on  the  target,  or  a reflec- 
tor mounted  on  the  target  to  return  light  generated  by  the 
tracker,  would  provide  a signal  to  the  star  tracker  whose 
outputs  would  indicate  the  angular  deviations  of  the  x^-axis 
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from  the  xLQS-axis. 


Pig.  14  Basic  STARTRAC  Concept  (Ref.  7«l-4) 

An  associated  problem  is  that  of  determining  the  atti- 
tude of  the  target  platform  relative  to  the  tracker’s  coord- 
inate frame.  A device  to  sense  the  angular  orientation  of 
the  platform  frame  with  respect  to  the  tracker  could  be  an 
autocollimator.  Such  a device  is  depicted  in  Fig.  15* 

In  an  autocollimator  a collimated  beam  of  light  is 
reflected  off  a mirror  on  the  target  and  the  return  beam  is 
passed  through  the  same  optics  and  focused  on  a photodetec- 
tor. The  displacement  at  the  detector  between  the  reflected 
image  and  the  effective  location  of  the  light  source  is  a 
measure  of  the  misalignment  between  the  units  (Ref.  9*17) • 
Harris,  et  al,  (Ref.  9*36-42)  discuss  two  systems  under 
development  to  accomplish  3-axis  optical  angular  measurements- 
the  Chrysler-Optical  Angular  Motion  Sensor  and  the  Barnes 


Flexure  Monitor  System.  Of  the  two  systems,  the  Chrysler 
system  is  reported  to  have  the  potential  of <5  sec  accuracy 
over  distances  greater  than  50  feet  and  the  Barnes  system  is 
said  to  have  accuracy  of  0.28  sec  over  a distance  of  5 feet 
(Ref.  907)*  A third  system  mentioned  briefly  in  Reference 
9 is  the  Martin  Universal  Star  Tracker  which  employs  a cor- 
relation tracker  using  stored  maps  of  the  light  sources  being 
tracked  (Ref.  9«40).  At  this  point,  the  tracker  hardware 
to  accomplish  range,  angle,  and  coordinate  frame  misalignment 
measurements  has  been  briefly  discussed.  A discussion  of  the 
accelerometers  in  the  system  follows. 


Fig.  15  Autocollimator  (Ref.  9*17) 

As  stated  by  Slater,  "The  integrating  type  of  accelero- 
meter is  particularly  attractive  for  gravity  measurements  on 
an  unstable  base  such  as  a ship  because  accelerations  associ- 
ated with  rolling  and  pitching  are  effectively  averaged  out." 
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(Ref.  24il96).  Slater  goes  on  to  point  out  that  bias  stabili- 
ty is  the  most  important  consideration  in  gravimetry.  A sub- 
stantial bias  component  in  accelerometer  output  would  effec- 
tively mask  out  any  gravitational  disturbance  variations 
unless  the  bias  could  be  compensated  for.  Although  in  this 
study,  accelerometer  errors  are  modelled  as  white  Gaussian 
noise,  it  is  recognized  that  any  operational  system  based  on 
this  proposal  would  probably  augment  the  state  vector  to 
include  estimating  the  bias, scale  factor  errors,  and  other 

non-linearities  in  the  accelerometers. 

<• 

In  addition  to  the  above  considerations,  it  should  be 
recognized  that  acceleration  components  associated  with  the 
vehicle’s  movement  through  the  air  or  water  must  be  known  to 
the  same  degree  of  accuracy  as  desired  for  the  gravity  dis- 
turbances. Slater  in  Reference  24  published  in  1964,  states 
that  accelerometer  sensitivities  of  10”  m/sec  had  been  dem- 
onstrated for  earth- based  gravimeters.  This  is  on  the  order 
of  the  anomalous  relative  gravitational  accelerations  encoun- 
tered in  this  study,  and  offers  encouragement  for  the  prac- 
tical application  of  the  gradient  mapping  system  proposed  here 

This  concludes  the  brief  discussion  of  the  hardware  which 
might  be  used  in  a practical  system  as  envisioned  here.  The 
EKF  formulation  is  discussed  next. 


III.  Extended  Kalman  Filter 


Formulation  (Ref.  15) 

It  is  assumed  that  the  reader  has  an  understanding  of 
linear  Kalman  filters  and  is  familiar  with  the  filter  equa- 
tions. Consequently,  the  discussion  which  follows  will  begin 
with  a statement  of  the  state  dynamics  model  for  the  EKF  and 
will  conclude  with  a description  of  the  EKF  equations.  No 
attempt  is  made  to  derive  the  EKF  formulation,  but  the 
interested  reader  is  referred  to  Ref.  6 or  Ref.  15  for  such 
details . 

The  EKF  is  one  method  used  in  estimation  problems  where 
nonlinearities  in  the  deterministic  portion  of  the  state 
dynamics  and  measurement  models  cannot  be  neglected.  In 
this  study,  the  state  equations  are  continuous  in  time,  but 
measurements  will  be  incorporated  at  discrete  intervals. 

For  this  reason,  the  continuous  form  of  the  EKF  with  discrete 
measurement  updates  is  appropriate. 

Let  the  system  state  equations  satisfy  the  following 
nonlinear  vector  stochastic  differential  equation  written  in 
the  familiar  white  noise  notation  as 


x(t)  = f +§(t)w(t) 


where i 

xfl. 


time  derivative  of  the  system  state 
n- vector  at  time  t, 


(3-D 


f[x(V] 


n- vector  of  nonlinear  system  dynamics. 


SCft-  system  noise  input  matrix,  and 
Y/W  = Gaussian  white  noise  vector  with  statistics 

E $ *(+')}  = g 

E £y;  W WT(t+'r)j  = Q(t)S  Cy) 

The  initial  conditions  for  the  state  vector  xC+.'i  is  assumed 

A 

to  be  a Gaussian  random  vector  with  mean  X0  and  covariance 
matrix  . The  carat  symbol,  or  "hat,"  over  a quantity  rep- 
resents an  estimate  of  that  quantity,  not,  in  general,  its 
true  value. 


Let  the  discrete- time  measurement  equations  satisfy  the 
following  nonlinear  vector  equation 


where  t 


= m- vector  of  measurements  at  time  "h  , 

] = vector  of  nonlinear  functions  of  states 

and  current  timeTi  , and 
= Gaussian  white  noise  vector  with  statistics 


With  the  system  dynamics  and  measurement  models  given  in  eqns 
(3-1)  and  (3-2),  one  desires  an  optimal  state  estimate.  In 
the  case  of  linear  systems  this  optimal  estimate  can  be  shown 
to  be  that  of  the  Kalman  filter.  To  extend  the  Kalman  filter 


are  X®  and  Po  ).  In  eqn.  (3-*0  above,  f [£(t/t,),t]  is 


defined  as  the  n-by-n  partial  derivative  matrix i 


F[S(t/tiU]  4 


i L-v'  '» J — ^ 

After  integrating  eqns.  (3-3)  and  (3-*0  to  time  'ft+i  , 

are  defined  asi 

8(ti«)  = £(Wh) 
f (t.vT)  = P (ti«/ti) 

and  used  in  the  next  measurement  update. 


(3-7) 


(3-8) 


(3-9) 


The  system  measurements, 2;  , are  incorporated  in  the  EKF 


k (ti)  = f (to  >8) 

5 (Til  -hl$(l0)igj  (3-9) 

p(ti+)=  fcrr) -K(^HLx(rr)jTi]p(Tr)  (3-10) 

where  tiKCTDjtn  is  defined  as  the  m-by-n  partial  deriva- 


(3-9) 

(3-10) 


tive  matrix i 


„[8Mir0a  Jfcgual^ 

and  K(T5  is  the  n-by-m  Kalman  gain  matrix.  With  the  basic 
EKF  formulation  completed,  one  can  turn  his  attention  to 
applying  this  non-linear  estimator  to  the  problem  under  study. 
Application  to  this  Problem 

The  state  equations  upon  which  the  estimator  is  based 
are  restated  here  using  slightly  different  notation  from 


that  used  in  Chapter  II.  The  9-state  vector,  X.(+)  , is  given 
by  (dropping  the  time  argument  for  convenience) 


where  X is  modelled  as  a Gaussian  random  vector  with  initial 
conditions  and  . The  non-linear  state  dynamics  can  be 
written  as  rr  T 


while  the  additive  state  dynamics  driving  noise  term 


Gffiwft.  is  given  by 


G(t)w(+)  = 


O 

w, 

W3 

Ws 

W* 

W7 


(3-14) 


l= 


where  G(t)  is  a 9x9  identity  matrix  and  wCt)  is  a 9-vector 
with  a zero  as  the  first  element.  The  w(t)  vector  has  sta- 
tistics 

E{w(f)}  = 0 
E^w(t)vv/T(t4-T^=  Q(t)  S(r) 


(3-15) 


where  a(f)  is 


Q(+)  = 


o 


V*. 

v»  o 

Vi  — 

<UST 

O S-Wp 

Vi 


<Lm 


(3-16) 


o 


and  SCT)  is  the  Dirac  delta  function.  Numerical  values  used 
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in  the  filter  for  X®  , jj*  , and  Q(+)  can  be  found  in  Appen- 
dix B. 

The  measurements,  as  discussed  in  Chapter  II,  are  linear 
measurements  of  states  Xt  • X5  , andXc,  and  satisfy  the 
equation 


2(t0  = H x(ti)  v(ti) 


(3-17) 


For  this  problem,  H x(ti)  is  expressed  as 


H x(t;)  = Xstri)  I (3-18) 

and  the  measurement  corruption  noise  vector  y(ti)  has  sta- 
tistics 

E£v(ti)}=0 

L 2 >Ti?=T} 

where  RCtO  is  given  by 


r,»  o o 
R(ti)  = | o ix*.  o 
O o 


(3-20) 


Numerical  values  for  R(t.)  are  also  given  in  Appendix  B. 

Finally,  two  additional  matrices  must  be  defined  before 


the  EKF  can  be  used  in  a simulation.  These  matrices  are 


BEWO»iQ  defined  in  eqn.  (3-11)  and  F[£(tA),+3  defined 
in  eqn.  (3-7).  For  this  problem  can  be  ex- 

pressed as  that  shown  in  Figure  16,  and 


H = 


\ oooooooo 

OOOOI oooo 

ooooo \ ooo 


(3-22) 


With  all  of  the  pieces  developed  that  are  needed  to 
operate  the  EKF,  the  simulation  begins  by  propagating  eqns. 
(3-3)  and  (3-4)  forward  one  time  interval  (10  3ec  was  used 
here).  The  10  sec  interval  was  chosen  because  the  data 
was  available  at  10  sec  intervals.  The  propagation  of  the 
state  estimate,  x(t)  , and  covariance  matrix,  f(t)  , was 
accomplished  by  direct  numerical  integration  of  eqns.  (3-3) 
and  (3-4).  The  integration  package  used  was  ODE,  an  ordinary 
differential  equation  solver  which  is  part  of  the  CDC  6600 
library  of  subroutines  at  the  Aeronautical  Systems  Division 
Computer  Center.  ODE  utilizes  a modified  divided  difference 
form  of  the  Adams  Pece  formulas  and  local  extrapolation.  It 
adjusts  the  order  and  step  size  to  control  the  local  error 
per  unit  step.  The  interested  reader  is  referred  to  Refer- 
ence 20  for  more  details  or.  this  computer  program. 

At  the  end  of  the  first  time  interval,  ODE  returns  the 
state  vector  and  covariance  matrix,  x(Tf)  and  fCtr)  respec- 
tively. At  this  point  system  measurements  are  taken.  These 
measurements  are  incorporated  into  the  EKF  and  used  to  improve 
the  state  and  covariance  estimates.  Equations  (3-8),  (3-9)* 
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and  (3-10)  are  typically  used  in  measurement  updating,  how- 
ever, for  this  problem  an  alternate  form  of  eqn.  (3-10)  was 
selected  initially.  This  alternate  form  of  the  covariance 
matrix  update  equation  is  sometimes  called  the  "Joseph  form" 
and  is  given  by  (Ref.  6«305) 

P't’-(l-Kd)  P'(l-  K Hf  + KRKT  (3-«) 

where  the  time  argument  is  again  dropped  for  convenience . 
Although  eqn.  (3-23)  involves  more  computation  time  than  eqn. 
(3-10),  it  is  inherently  more  accurate  and  it  was  on  this 
basis  that  the  Joseph  form  was  selected  initially  for  this 
filter.  Unfortunately,  despite  its  accuracy,  covariance 
measurement  updating  by  eqn.  (3-23)  does  not  insure  retaining 
positive  semidefiniteness  of  the  error  covariance.  This  loss 
of  positive  semidefiniteness  was  observed  during  the  first 
few  simulation  runs.  As  stated  by  Maybeck  (Ref.  I6tl2),  it 
can  be  shown  that  the  Kalman  filter  algorithm  is  numerically 
unstable.  Truncation  errors  due  to  finite  wordlength  are 
thus  propagated  once  they  occur  and  do  not  die  out.  In  this 
case,  a CDC  6600  computer  with  a 6o-bit  wordlength  was  usedt 
however,  even  with  this  large  a wordlength  numerical  insta- 
bility in  the  algorithm  was  observed.  Such  problems  often 
occur  when,  as  Maybeck  puts  iti 

"...  ( 1 ) the  measurements  are  very  accurate 
(eigenvalues  of  gCtt)  are  small  relative  to 
those  of  f(tr)  , this  being  accentuated  by 
large  eigenvalues  in JE*  ) or  (2)  a linear 
combination  of  state  vector  components  is 
known  with  great  precision  while  other 
combinations  are  nearly  unobservable  (i.e., 
there  is  a large  range  of  magnitudes  of 
state  covariance  eigenvalues."  (Ref.  16i12,13) 
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Certainly  in  this  application  condition  (1)  above  applies. 
What,  then  are  the  possible  solutions  to  this  problem  ? 

Traditional  solutions  tc  the  inherent  numeric  problems 
of  the  Kalman  filter  include  performing  measurement  updating 
in  double  precision,  scaling  variables  to  reduce  dynamic 
range,  or  employing  ad  hoc  modifications  (e.g.,  placing  lower 
limits  on  covariance  diagonal  terms  or  upper  limits  on  off- 
diagonal  terms).  In  this  author's  mind,  a solution  prefer- 
able to  those  just  mentioned  is  to  employ  one  of  the  various 
"square  root  filters"  which  have  been  developed. 

The  basic  idea  behind  square  root  filtering  is  to  prop- 
agate and  update. an  error  covariance  square  root  or  inverse 
covariance  square  root  rather  than  the  covariance  or  its 
inverse  themselves.  The  advantage  of  this  approach,  again 
quoting  Maybeck,  include  the  following 

"The  square  root  approach  can  yield  twice  the 
effective  precision  of  the  conventional 
filter  in  ill-conditioned  problems.  In 
other  words,  the  same  precision  can  be 
achieved  with  approximately  half  the  wordlength. 
Moreover,  this  method  is  completely  success- 
ful in  maintaining  the  positive  semidefinite- 
ness of  the  error  covariance."  (Ref.  I61I3) 

Among  the  square  root  recursions  that  have  been  developed 
are  the  Potter  covariance  square  root,  Carlson  covariance 
square  root,  inverse  covariance  square  root,  and  the  U-D 
covariance  factorization.  The  reader  is  referred  to  Reference 
16  for  an  excellent  discussion  of  these  four  square  root 
recursive  algorithms. 

The  approach  finally  adopted  for  this  problem  was  the 
U-D  covariance  factorization  method.  Although  not  strictly 


a square  root  filter,  the  U-D  covariance  factorization  filter 
is  closely  related  and  provides  all  of  the  numerical  pre- 
cision of  the  square  root  filters  but  without  performing  any 
time  consuming  square  root  operations.  Actually,  for  this 
problem  the  U-D  factorization  was  only  used  for  the  covari- 
ance and  state  estimate  updating.  Propagation  of  the  state 
estimate  and  covariance  between  measurements  was  performed 
by  numerically  integrating  eqns.  (3-3)  and  (3-4)  as  discussed 
previously.  This  was  done  because  it  was  felt  that  only 
during  measurement  updating  was  the  numerical  instability  of 
the  Kalman  filter  algorithm  evidenced,  and  because  U-D  prop- 
agations depend  on  the  linearity  of  the  model.  The  U-D 
covariance  factorization  recursion  for  updating  the  state 
estimate  and  error  covariance  will  bow  be  given.  No  attempt 
is  made  to  derive  the  equations  but  only  to  give  the  filter 
algorithm  itself.  For  additional  details  and  discussion  see 
Reference  16,  section  8.  What  follows  was  taken  from  this 
same  source. 

U-D  Covariance  Factorization  Update  (Ref.  16 » 42-47) 

This  method  expresses  the  covariance  before  and  after 
measurement  updating  as 


P ' = U'D'  U-T 

(3-24) 

£*  = U+D+U+T 

(3-25) 

where  U matrices  are  upper  triangular  and  unitary  (i.e., 
ones  along  the  diagonal)  and  the  D matrices  are  diagonal. 
As  was  stated  previously,  the  U-D  filter  is  not  strictly 
speaking  a square  root  filter,  but  it  is  considered  such 
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because  UD*  corresponds  directly  to  the  covariance  square 
root. 

The  algorithm  as  employed  here  begins  at  the  end  of  the 
current  time  propagation  interval  * when  P(tr)  , a nxn 

symmetric,  positive  semidefinite  matrix  is  returned  by  ODE. 
The  first  step,  then,  is  to  generate  the  U and  D factors 
of  JP  . The  following  algorithm  will  generate  a uniquely 
defined  pair  of  U-D  factors. 

First,  for  the  n-th  column 


Dnn  fn 


Uin 


te-n 


(3-26a) 


P*n/Dnn  ia  jl 


Then  for  the  remaining  columns,  for  j = n-1,  n-2, 


compute i 


Pjj  = *3l3  “2  Dkk  Uj* 

K=J+l 

Uij"[cn.-iDKKUi»uj, 

A A • _ 


i»j 


v%y  — 1 • n / ' “*  (3-26b) 

( [Rj  Dkx  U i*,  UjkJ/Du 

The  matrices  U~and  £f have  now  been  computed.  The  scalar 
measurement  update  for  the  U-D  filter  will  now  be  given. 

Using  the  measurement  value  and  the  known  I XH 
and  scalar  f(Tf)  , one  computes 

f - utth.t  ~) 


vj = Pjjfs 

cu  = r 


j — 


(3-27) 


Then  for  K « 1,  2,  ...,  n,  one  calculates 
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Qk  — <Xk.-|  + 'fiiVK 

= &KK  aK-»/aK 

Pk  =~^K/aK-l 
UjK  = HjV  + bj  Pk 
bj  **-  bj  + IAjk  vk  1 


»K  =UjK+bjPK  ■) 

bj  — bj  + U"kvk  j j->A~»(*-0 


(3-28) 


In  eqn.  (3-28),  "^J*  denotes  replacement,  or  "writing  over" 

old  variables  for  efficiency.  For  K=l,  only  the  first  three 
equations  are  needed.  After  the  n iterations  of  eqn.  (3-28), 
U*  and  D*  have  been  computed,  and  the  filter  gain  K(-h) 
can  be  calculated  in  terms  of  the  n- vector  b made  up  of  com- 
ponents b|  * bi . •«.,br\  computed  in  the  last  iteration  of 
eqn.  (3-28),  and  the  state  updated  as 

K(+V)  = b /an 

Vector  measurements  updates  for  RCtf)  diagonal,  as  was  the 
case  here,  are  performed  component  by  component. 

The  covariance  update  is  computed  by  eqn.  (3-25)  and 
then  ODE  is  called  again  to  propagate  the  state  and  covari- 
ance ahead  to  the  next  measurement  time  (ti.O  . 

Before  leaving  this  subject,  a brief  comparison  owed  to 
Maybeck  (Ref.  16« 5^-57)  will  be  given.  Table  III  presents 
the  number  of  operations  required  for  one  time  propagation 
and  one  measurement  update  for  the  case  of  n*10,  s=10,  and 
ma2,  where  n is  the  state  dimension,  s is  the  noise  dimension. 


8 is  the  noise  dimension,  and  m is  the  measurement  dimension 


The  last  column  of  Table  III  gives  the  processor  time  for  one 
filter  cycle  based  on  the  following  instruction  times  typical 
of  the  IBM  360 i 


addition 


multiplication 

division 


square  root 


Table  III.  Operations  for  One  Filter  Cycle 

(n=s=lO,  m=2) 


SQUARE 


TIME 


ADDS 


MULTIPLIES 


DIVIDES 


Conventional 

Kalman 


Joseph  Form 
Kalman 


As  Table  III  shows,  for  the  case  given,  the  U-D  filter 
falls  between  the  conventional  Kalman  filter  and  the  Joseph 
form  Kalman  filter  in  processing  time.  However,  for  a modest 
increase  in  computational  burden  over  the  conventional  Kalman 
filter  one  avoids  its  numerical  instability  and  gains  twice 
the  accuracy  in  the  same  wordlength  by  employing  the  U-D 
filter.  Based  on  Table  III,  one  wonders  why  it  would  ever  be 
advantagous  to  use  the  Joseph  form  over,  say,  the  U-D  covari- 
ance factorization. 

The  initial  choice  of  the  Joseph  form  for  use  in  the 


problem  under  study  here  was  motivated  by  a desire  for 
numerical  accuracy.  What  was  overlooked  was  the  fact  that 
numerical  accuracy  does  not  insure  numerical  stability. 

Truth  Model 

The  truth  model  for  this  problem  consisted  of  a 4-state 
model  of  the  system  dynamics  plus  the  3 deterministic  gravity 
disturbance  states.  The  true  state  vector  is  given  by, 


and  the  state  differential  equations  are 


x»  (xl+xj)  + S3 

(-2-XxXs  -Sg^/Xi 
(-iX^X.  +S3yVXi 


In  the  truth  model,  perfect  tracking  was  assumed,  and  thus, 
Sn  and  &G  are  zero  for  all  time.  This  assumption  was  made 
since  the  tracking  loop  would  typically  be  operating  with  a 
much  lower  period  than  the  filter  cycle  time  of  10  sec.  In 
view  of  this  it  was  decided  not  to  model  tracker  dynamic?. 

In  any  further  work  in  this  area,  a shorter  filter  cycle 
could  be  used  and  a dynamic  model  of  the  tracker  added  to  the 


truth  model. 

The  ^ components  were  determined  by  runnir^  the  program 
GTGRY  discussed  in  Reference  25*  Inputs  to  GTGRV  were  data 
cards  with  the  positions  (expressed  in  e-frame  coordinates) 
at  which  the  gravity  vector  ^ and  gradient  matrix  _P  were 
desired.  Formulas  used  by  GTGRV  to  calculate  the  gravity 
vector  and  gradient  contributions  of  each  point  mass  in  the 
gravity  field  are  given  in  Reference  25,  Appendix  B,  pp.  5-7. 
Seven  point-mass  sets  were  used  in  the  calculations  with  grid 
sizes  of  5°,  1°,  15’,  5',  2.5’,  1.25’,  and  0.3125',  where 
"(o)"  means  degree  of  longitude  and  latitude  and  means 

minute  of  longitude  and  latitude  (60'=1°).  The  outputs  of 
GTGRV  used  in  this  study  were  the  % and  Pxx  , Pxy  , and 
Fxe  components  expressed  in  e-frame  coordinates. 

The  "true"  state  vector,  Xt  , was  obtained  by  direct 
numerical  integration  of  the  differential  equations  in  eqn. 
(3-31)  using  the  % components  calculated  by  GTGRV.  There 
were  no  noise  inputs  to  the  state  dynamics. 

Gravity  Gradient  Estimate 

The  gravity  gradient  JH  as  discussed  in  Appendix  A can 


be  written  as 


r = ii_ 


(3-32) 


where, is  the  gravity  vector  as  a function  of  position, 
JT  , and  to  first  order, 


■|ltU|kl=  fiL  II 

_5rJ  &£  a Sry 


(2-10) 


Obviously,  the  accuracy  with  which  one  can  approximate 
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the  gradient  with  this  technique  is  directly  proportional  to 
the  accuracy  of  the  ZQ.  estimate.  Likewise,  one  would  expect 
a better  approximation  to  the  true  P for  Id  small . Of 
course,  as  lei  decreases  so  does  |ss(r)|  , which  poses  prob- 
for  the  EKF  estimator  in  view  of  the  noise  process  strengths. 
Figures  17  and  18  show  plots  of  the  true  fix  * f"*y » 
and  Pxz  components  for  both  altitude  cases  simulated. 

Notice  that  the  true  P values  are  in  the  range  of  several 
tens  of  EU  (Eotvo's  Units  = 10~^sec"2) . In  this  simulation, 
only  a tracker- target  pair  along  the  x -axis  was  used.  Thus, 

V 

the  components  estimated  by  eqn.  (2-10)  for  this  problem  are 
fxx*  Pxy  * and  Px^  , or  the  first  column  of  the_P  matrix. 
With  one  pair  of  trackers  and  targets  one  gets  three  of  the 
five  independent  elements  of  P . One  additional  pair  ortho- 
gonal to  the  first  pair,  say,  along  the  y -axis  would  provide 

c 

the  other  two  independent  elements  and  the  entire  gravity 
gradient  tensor  would  be  determined.  In  practice,  a third 
pair  along  the  z -axis  would  probably  be  desirable  so  that 
data  smoothing  could  be  exploited  to  hopefully  improve  the 
estimate  of£  . 

The  added  components  in  the  full-up  system  instrumenting 
all  three  axes  would  not  increase  the  state  vector  dimension 
by  a factor  of  three.  Instead,  the  additional  two  axes  would 
add  two  range  states,  two  range  rate  states,  one  state 

and  one  angular  deviation  state. 


Monte  Carlo  Analysis  and  Results 


Monte  Carlo  Method 


The  Monte  Carlo  study  is  one  of  two  fundamental  tech- 
niques of  conducting  a performance  (sensitivity)  analysis  of 
a prospective  Kalman  filter.  One  seeks  an  accurate  statisti- 
cal portrayal  of  the  filter’s  estimation  errors  without 
having  to  build  the  filter  and  test  it  in  a "real  world" 
environment.  One  can  do  this  via  computer  simulation  where- 
in the  "real  world"  is  replaced  by  a truth  model,  the  most 
complete  mathematical  model  one  can  devise  to  simulate  the 
real  world. 

Figure  19  illustrates  how  a performance  analysis  of  a 
Kalman  Filter  might  be  conducted.  This  figure  does  not  rep- 
resent the  most  general  performance  evaluation  mechanization 
but  is  tailored  to  the  needs  of  the  problem  under  study. 


Kalman 

Filter 


Fig.  19  Performance  Evaluation  of  a Kalman  Filter 


Prom  the  truth  model,  the  true  state  vector,  *r  • and 
the  true  covariance, £r  , are  generated.  The  white  Gaussian 
noise  process,  >£*  , drives  the  truth  model  to  accurately 
generate  the  measurements,  Z . The  Kalman  filter  combines 
the  measurement  information  with  its  own  internal  dynamics 
model  of  the  system  and  produces  its  own  "best"  estimate  of 
the  states,  &,  and  covariance,  F . The  error  process, » 
and  error  covariance,  £*  , are  formed  by  the  "difference" 
between  the  true  states  and  the  filter  states.  This  "differ- 
ence" must  be  modified  from  its  accepted  meaning  to  account 
for  the  differences  in  dimension  between  the  truth  model  and 
filter  model  in  the  general  case.  The  objective  of  the  per- 
formance analysis  is  to  statistically  characterize,  §T  . 

Already  mentioned  as  a technique  for  generating  the 
statistics  of  is  the  Monte  Carlo  study.  The  other  widely 
used  technique  is  a mean  and  covariance  analysis.  This 
technique  is  exploited  when  the  truth  model  is  a linear  sys- 
tem driven  by  white  Gaussian  noise  and  having  linear  measure- 
ments corrupted  by  white  Gaussian  noise.  Since  a non-linear 
system  is  involved  in  this  study,  the  covariance  analysis 
is  not  used  and  will  not  be  discussed  further. 

In  conducting  a Monte  Carlo  study,  many  samples  of  gT 
are  generated  by  simulation,  and  the  ensemble  average,  or 
statistical  behavior  of  &T  is  computed  directly.  Theoretica- 
lly, an  infinite  number  of  simulations  is  necessary  to  insure 
that  the  computed  sample  statistics  do,  in  fact,  portray  the 
filter's  performance.  Obviously,  conducting  an  infinite 
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number  of  simulations  is  out  of  the  question,  so  only  a 
limited  number  of  simulations  are  conducted.  Choosing  the 
number  of  simulations  which  will  result  in  process  statistics 
which  adequately  reflect  the  "true"  statistics  is  often  in- 
fluenced by  the  need  to  conserve  computer  time  and  stay  with- 
in budget  constraints.  Although  there  are  no  rules  which 
govern  the  choosing  of  am  appropriate  number  of  simulations, 
the  degree  of  non-linearity  in  the  truth  model  is  certainly 
a factor.  For  problems  with  insignificant  non-linearities 
and  modest  state  dynamics  perhaps  as  few  as  5-10  simulations 
would  suffice.  On  the  other  hand,  for  highly  non-linear 
models  as  many  as  15-25  simulations  (or  more)  may  be  desir- 
able. Again,  time  and  budget  constraints  may  well  dictate 
an  upper  limit,  but  one  should  be  wary  of  drawing  conclusions 
about  filter  performance  based  on  inadequate  statistics. 

One  means  of  choosing  an  appropriate  number  of  simula- 
tions is  to  plot  each  states'  changes  in  computed  variance 
versus  the  number  of  simulations.  Such  a plot  will  typically 
reach  nearly  steady  state  conditions  where  adding  more  simu- 
lations does  little  to  the  sample  statistics.  An  appropriate 
number  of  simulations  can  then  be  chosen  from  the  steady 
state  region. 

The  Monte  Carlo  analysis  was  never  completed  in  this 
case  because  of  severe  numerical  problems  encountered  during 
the  simulations.  These  problems  are  discussed  next. 

Results 

The  Monte  Carlo  analysis  was  never  successfully 
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completed i and,  consequently,  few  results  can  be  reported. 

Prom  the  beginning  of  the  start  of  simulation  runs,  the  EKF 
exhibited  divergence.  Within  two  or  three  filter  cycles,  the 
P matrix  would  lose  positive  semidefiniteness  and  the 
state  estimates  would  begin  to  grow, seemingly  without  bound. 
Many  attempts  were  made  to  correct  the  problem. 

The  first  attempt  was  to  try  to  tune  the  filter.  Tuning 
a Kalman  Filter  involves  adjusting  the  _P«  , ClCf)  • and  RCtV) 
values  to  obtain  the  "best"  performance  possible.  Many 
tuning  runs  were  made  but  to  no  avail. 

This  problem  was  known  to  be  ill-conditioned  numerically 
because  of  the  wide  range  of  values  between  state  variables 
(range  was  on  the  order  of  10m  while  Sj.  components  got  as 
small  as  10  m/sec  ) . The  decision  to  change  the  form  of  the 
filter  update  equations  was  based  on  a desire  to  exploit 
the  superior  numerical  characteristics  of  the  U-D  covariance 
factorization  update  algorithm.  When  this  change  was  made, 
loss  of  positive  semidefiniteness  of  P was  still  observed 
although  it  did  not  occur  until  the  third  or  fourth  filter 
cycle . 

Once  again  many  tuning  runs  were  made  but  without  allevi- 
ating the  problem.  It  was  assumed,  at  this  point,  that  the 
numerical  ill-conditioning  of  the  problem  was  at  the  heart  of 
the  filter  instability. 

With  time  running  out,  a final  attempt  was  made  to  scale 
the  variables.  This  was  done  by  scaling  each  variable  to 
fall  in  the  range  of  0<X^2. 


li 
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The  formula  used  for  scaling  was 
v,  X — Xmin 


where  t 


maximum  value  of  X based  on  the 


truth  model  data,  and 


minimum  value  of  X based  on  the 


truth  model  data 


After  scaling,  several  tuning  runs  were  conducted.  It 
was  discovered  that  only  by  increasing  the  values  for  the 
elements  of  R(Ti)  (i.e.,  increasing  the  uncertainty  in  the 
measurements)  could  positive  semidefiniteness  in  P be  main- 
tained. The  original  RCti)  was  given  by 


while  the 


for  the  only  case  where  the  filter  worked 


successfully  was 


(4-3) 


The  inability  of  the  filter  to  avoid  numerical  instabili 
ty  even  after  scaling  and  the  U-D  covariance  update  were 
implemented  was  entirely  unexpected.  The  explanation  lies  in 
the  small  eigenvalues  of  RCn')  as  originally  proposed.  Mea- 
surements of  this  accuracy  appear  to  the  filter  to  be  nearly 
perfect,  and  this  leads  to  numerical  instability  in  the  Kal- 
man filter  algorithm.  This  indicates  that  perhaps  the  EKF 
was  a poor  choice  of  estimator  to  employ  in  this  problem  and 
leads  to  a number  of  recommendations  to  be  covered  in  Chapter 


V.  Conclusions  and  Recommendations 


Conclusions 


As  has  already  been  stated,  the  lack  of  positive  results 
in  this  study,  although  disappointing  in  the  extreme,  does 
not  invalidate  the  proposal  advanced  in  this  thesis.  Rather, 
it  has  been  shown  that  the  problem  is  extremely  ill-condi- 
tioned and  the  Extended  Kalman  Filter  is  not  an  appropriate 
choice  for  this  application.  As  a result  of  the  effort 
expended  in  the  course  of  this  study,  the  following  have  been 
achieved i 

(1)  A proposal  has  been  made  for  am 
alternative  approach  to  gravity 
gradient  estimation. 

(2)  The  state  equations  upon  which  to 
base  am  estimator  have  been  developed. 

(3)  This  study  delineated  some  of  the 
numerical,  problems  which  amy  sub- 
sequent work  in  this  area  must 
address.  Some  tentative  solutions 
were  explored  without  significant 
improvement  in  the  results. 

(4)  The  Extended  Kalman  Filter  has 
been  shown  to  be  a poor  choice 
for  this  application. 

(5)  The  method  of  gravity  gradient 
estimation  proposed  here  should 
probably  be  restricted  to  use 

in  a ship  or  submarine  environment. 

Recommendations 

First,  a plea  must  be  made  for  someone  to  continue  the 
investigation  started  here.  The  proposed.  deserves  am  adequate 
test  of  its  feasibility. 

The  numerical  problems  encountered  in  this  study  should 
be  a lesson  for  anyone  attempting  a similar  study.  The 
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the  gradient  estimate. 


Finally,  the  extension  of  this  study  to  a full  three 
axis  implementation  would  be  valuable  in  ascertaining  the 
practical  limits  to  the  gradient  estimation  accuracies  that 
could  be  achieved  by  this  technique.  The  potential  future 
benefits  of  more  research  in  this  area  include  better  know- 
ledge  of  the  earth's  gravity  field  and  subsequent  improve- 
ments in  aerospace  and  seaborne  vehicle  navigation  and 
guidance . 


o 
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Appendix  A 


The  Gravity  Gradient  Tensor 


A great  deal  of  research  is  being  conducted  today  on  a 
new  class  of  inertial  instruments,  gradiometers , whose  purpose 
is  to  measure  elements  of  the  gravity  gradient  matrix.  The 
system  studied  elsewhere  in  this  report  is  another  related 
approach  to  accomplish  the  same  objective.  What  follows  is 
a brief  discussion  of  the  mathematical  properties  of  the 
gravity  gradient  tensor. 

It  is  well  known  that  gravitational  potential  is  a scalar 
function  <f>  of  position  H relative  to  the  attracting  body 
(see  for  example  Reference  10).  The  gravitational  acceleration 
is  the  spatial  derivative,  or  gradient,  of  the  gravita- 
tional potential  function.  Expressed  in  Cartesian  coordinates, 


The  gravity  gradient  tensor 
derivative  of  . Expre; 


then 


Since  the  gravity  field  is  conservative,  we  have  from 


Equation  A- 3 expresses  gravity  gradient  tensor  symmetry  and 
equation  A-4, known  as  Laplace's  equation, requires  the  trace 
of  the  gradient  matrix  to  be  equal  to  zero.  These  two 
considerations  reduce  the  number  of  independent  elements  of 
the  gradient  matrix  to  five.  Thus,  complete  determination 
of  the  tensor  at  any  point  requires  five  independent  measure 
ments . 

In  keeping  with  convention,  the  symbol  j]1  is  used  to 
represent  the  gravity  gradient  tensor, 


Appendix  B 

Numerical  Values  for  , P*  , , and  RCrO 

The  following  are  the  numerical  values  employed  in 
the  EKF  used  in  this  study. 
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Appendix  C 

PSD  Plots  for  Disturbance  Data 
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The  following  six  plots  were  obtained  from  the  150 
data  points  representing  the  true  gravity  differential  dis- 
turbance data.  The  IMSL  library  subroutine  FTFREQ  was  used 
to  compute  the  power  spectra  for  plotting  on  the  CALCOMP 
plotter. 
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PSD  Plot  of  Delta  G 
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